home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / LIBRARY / SHDK_2 / SHCMPLX.PAS < prev    next >
Pascal/Delphi Source File  |  1992-05-13  |  13KB  |  485 lines
  1. {$I SHDEFINE.INC}
  2.  
  3. {$I SHUNITSW.INC}
  4.  
  5. {$D-,L-}
  6. unit ShCmplx;
  7. {
  8.                                 ShCmplx
  9.  
  10.                        A Complex Arithmetic Unit
  11.  
  12.                                    by
  13.  
  14.                               Bill Madison
  15.  
  16.                    W. G. Madison and Associates, Ltd.
  17.                           13819 Shavano Downs
  18.                             P.O. Box 780956
  19.                        San Antonio, TX 78278-0956
  20.                              (512)492-2777
  21.                              CIS 73240,342
  22.  
  23.                   Copyright 1991 Madison & Associates
  24.                           All Rights Reserved
  25.  
  26.         This file may  be used and distributed  only in accord-
  27.         ance with the provisions described on the title page of
  28.                   the accompanying documentation file
  29.                               SKYHAWK.DOC
  30. }
  31.  
  32. interface
  33. uses
  34. TPCRT,
  35. TPDOS,
  36. TPSTRING,
  37.   TpInline,
  38.   ShErrMsg;
  39.  
  40. type
  41. {$IFNDEF Gen87}
  42.   extended        = real;
  43. {$ENDIF}
  44.   ComplexElement  = extended;
  45.   ComplexBaseType = record
  46.                       Re,
  47.                       Im  : ComplexElement;
  48.                       end;
  49.   Complex         = ^ComplexBaseType;
  50.  
  51.   PFreeRec        = ^TFreeRec;
  52.   TFreeRec        = record
  53.                       Next  : PFreeRec;
  54.                       Size  : pointer;
  55.                       end;
  56.  
  57. const
  58.   cxHaltErr     : boolean        = true;    {Stop program execution on
  59.                                               non-I/O errors}
  60.   cxOK          = 0;
  61.   cxDivByZero   = 200;
  62.   cxBadMagnitude= 201;
  63.  
  64. procedure CmplxDispose(A : Complex);
  65. {Replacement for the normal TP Dispose procedure. MUST be used to dispose
  66.  of variables of type Complex.}
  67.  
  68. procedure CmplxInit;
  69. {Initializes the complex arithmetic package.}
  70.  
  71. procedure CmplxDeInit;
  72. {De-initializes the complex arithmetic package, releasing all the ring
  73.  buffer heap space.}
  74.  
  75. function CmplxError  : integer;
  76.  
  77. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  78. {Converts a complex value to a string of the form (Re + Im i)}
  79.  
  80. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  81. {Converts a complex in polar form to a string with the angle in radians.}
  82.  
  83. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  84. {Converts a complex in polar form to a string with the angle in degrees.}
  85.  
  86. procedure CAbs(A  : Complex; var Result : ComplexElement);
  87. function CAbsF(A  : Complex)  : ComplexElement;
  88. {Returns the absolute value of a complex number}
  89.  
  90. procedure CConj(A : Complex; Result : Complex);
  91. function CConjF(A : Complex) : Complex;
  92. {Returns the complex conjugate of a complex number}
  93.  
  94. procedure CAdd(A, B : Complex; Result : Complex);
  95. function CAddF(A, B : Complex) : Complex;
  96. {Returns the sum of two complex numbers A + B}
  97.  
  98. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  99. function RxCF(A : ComplexElement; B : Complex) : Complex;
  100. {Returns the complex product of a real and a complex}
  101.  
  102. procedure CSub(A, B : Complex; Result : Complex);
  103. function CSubF(A, B : Complex) : Complex;
  104. {Returns the difference between two complex numbers A - B}
  105.  
  106. procedure CMul(A, B : Complex; Result : Complex);
  107. function CMulF(A, B : Complex) : Complex;
  108. {Returns the product of two complex numbers A * B}
  109.  
  110. procedure CDiv(A, B : Complex; Result : Complex);
  111. function CDivF(A, B : Complex) : Complex;
  112. {Returns the quotient of two complex numbers A / B}
  113.  
  114. procedure C2P(A : Complex; Result : Complex);
  115. function C2PF(A : Complex) : Complex;
  116. {Transforms a complex in cartesian form into polar form}
  117. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  118.  
  119. procedure P2C(A : Complex; Result : Complex);
  120. function P2CF(A : Complex) : Complex;
  121. {Transforms a complex in polar form into cartesian form}
  122. {The magnitude is stored in A^.Re and the angle in A^.Im}
  123.  
  124. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  125. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  126. {Raises a complex (in polar form) to a real power}
  127.  
  128. {===========================================================}
  129.  
  130. implementation
  131.  
  132. VAR OUTPUT : TEXT;
  133.  
  134. const
  135.   MaxStack            = 5;
  136.   StackTop            = MaxStack - 1;
  137.   Pi                  = 3.14159265358979;
  138.   PiOver2             = Pi / 2.0;
  139.   TwoPi               = Pi * 2.0;
  140.   F180OverPi          = 180.0 / Pi;
  141.   SpConj    : byte    = 0;
  142.   SpSum     : byte    = 0;
  143.   SpDiff    : byte    = 0;
  144.   SpRProd   : byte    = 0;
  145.   SpProd    : byte    = 0;
  146.   SpQuot    : byte    = 0;
  147.   SpPolar   : byte    = 0;
  148.   SpCartsn  : byte    = 0;
  149.   SpPower   : byte    = 0;
  150.   CmpInited : boolean = false;
  151.  
  152. type
  153.   StackArray  = array[0..StackTop] of Complex;
  154.  
  155. var
  156.   Conj,
  157.   Sum,
  158.   Diff,
  159.   RProd,
  160.   Prod,
  161.   Quot,
  162.   Polar,
  163.   Cartsn,
  164.   Power    : StackArray;
  165.   CmpError : integer;
  166.  
  167.  
  168. procedure CmplxDispose(A : Complex);
  169. {Replacement for the normal TP Dispose procedure. MUST be used to dispose
  170.  of variables of type Complex.}
  171.   var
  172.     HP      : pointer;
  173.     B       : PFreeRec;
  174.     IsFree  : boolean;
  175.     FRec    : TFreeRec;
  176.   begin
  177.     HP := Normalized(HeapPtr);
  178.     A := Normalized(A);
  179.     B := Normalized(FreeList);
  180.     while (A <> Complex(B)) and (B <> HP) do begin
  181.       B := Normalized(B^.Next);
  182.       end;
  183.     if (A = Complex(B)) then begin
  184.       exit;
  185.       end;
  186.     Dispose(A);
  187.     end; {CmplxDispose}
  188.  
  189. procedure ChkInit;
  190.   begin
  191.     if not CmpInited then
  192.       RunErrorMsg(400, 'Unit ShCmplx not initialized.');
  193.     end;
  194.  
  195. procedure CmplxInit;
  196.   var
  197.     T1  : byte;
  198.   begin {CmplxInit}
  199.     for T1 := 0 to StackTop do begin
  200.       New(Cartsn[T1]);
  201.       New(Conj[T1]);
  202.       New(Diff[T1]);
  203.       New(Polar[T1]);
  204.       New(Power[T1]);
  205.       New(Prod[T1]);
  206.       New(Quot[T1]);
  207.       New(RProd[T1]);
  208.       New(Sum[T1]);
  209.       end;
  210.     CmpInited := true;
  211.     end; {CmplxInit}
  212.  
  213. procedure CmplxDeInit;
  214. {De-initializes the complex arithmetic package, releasing all the ring
  215.  buffer heap space.}
  216.   var
  217.     T1  : byte;
  218.   begin {CmplxDeInit}
  219.     for T1 := StackTop downto 0 do begin
  220.       {Flush the heap space}
  221.       if    Sum[T1] <> nil then
  222.         Dispose(Sum[T1]);
  223.       if  RProd[T1] <> nil then
  224.         Dispose(RProd[T1]);
  225.       if   Quot[T1] <> nil then
  226.         Dispose(Quot[T1]);
  227.       if   Prod[T1] <> nil then
  228.         Dispose(Prod[T1]);
  229.       if  Power[T1] <> nil then
  230.         Dispose(Power[T1]);
  231.       if  Polar[T1] <> nil then
  232.         Dispose(Polar[T1]);
  233.       if   Diff[T1] <> nil then
  234.         Dispose(Diff[T1]);
  235.       if   Conj[T1] <> nil then
  236.         Dispose(Conj[T1]);
  237.       if Cartsn[T1] <> nil then
  238.         Dispose(Cartsn[T1]);
  239.       {Reset the heap pointers}
  240.       Conj[T1] := nil;  Sum[T1] := nil;     Diff[T1] := nil;
  241.       RProd[T1] := nil; Prod[T1] := nil;    Quot[T1] := nil;
  242.       Polar[T1] := nil; Cartsn[T1] := nil;  Power[T1] := nil;
  243.       end;
  244.     {Reset the ring buffer pointers}
  245.     SpConj := 0;  SpSum := 0;     SpDiff := 0;
  246.     SpRProd := 0; SpProd := 0;    SpQuot := 0;
  247.     SpPolar := 0; SpCartsn := 0;  SpPower := 0;
  248.     CmpInited := false;
  249.     end; {CmplxDeInit}
  250.  
  251. function CmplxError  : integer;
  252.   begin
  253.     CmplxError := CmpError;
  254.     CmpError := cxOK;
  255.     end;
  256.  
  257. function Real2Str(A : ComplexElement; Width, Places : byte) : string;
  258.   var
  259.     T1  : string;
  260.   begin
  261.     Str(A:Width:Places, T1);
  262.     Real2Str := T1;
  263.     end;
  264.  
  265. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  266.   begin
  267.     if A^.Im >= 0.0 then
  268.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
  269.                        Real2Str(A^.Im, Width, Places) + 'i)'
  270.     else
  271.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
  272.                        Real2Str(Abs(A^.Im), Width, Places) + 'i)';
  273.     end;
  274.  
  275. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  276. {Converts a complex in polar form to a string with the angle in degrees.}
  277.   begin
  278.     CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
  279.       Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
  280.     end;
  281.  
  282. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  283. {Converts a complex in polar form to a string with the angle in radians.}
  284.   begin
  285.     CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
  286.       Real2Str((A^.Im),Width,Places)+' rad)';
  287.     end;
  288.  
  289. procedure IncPtr(var StackName : StackArray; var StackPtr : byte);
  290.   begin
  291.     StackPtr := (StackPtr + 1) mod StackTop;
  292.     if StackName[StackPtr] <> nil then begin
  293.       Dispose(StackName[StackPtr]);
  294.       end;
  295.     New(StackName[StackPtr]);
  296.     end;
  297.  
  298. function CAbsF(A  : Complex) : ComplexElement;
  299.   begin
  300.     CmpError := cxOK;
  301.     CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
  302.     end;
  303.  
  304. procedure CAbs(A  : Complex; var Result : ComplexElement);
  305.   begin
  306.     Result := CAbsF(A);
  307.     end;
  308.  
  309. function CConjF(A : Complex) : Complex;
  310.   begin
  311.     CmpError := cxOK;
  312.     ChkInit;
  313.     Conj[SpConj]^.Re := A^.Re;
  314.     Conj[SpConj]^.Im := -A^.Im;
  315.     CConjF := Conj[SpConj];
  316.     IncPtr(Conj, SpConj);
  317.     end;
  318.  
  319. procedure CConj(A : Complex; Result : Complex);
  320.   begin
  321.     Result^ := CConjF(A)^;
  322.     end;
  323.  
  324. function CAddF(A, B : Complex) : Complex;
  325.   begin
  326.     CmpError := cxOK;
  327.     ChkInit;
  328.     ChkInit;
  329.     Sum[SpSum]^.Re := A^.Re + B^.Re;
  330.     Sum[SpSum]^.Im := A^.Im + B^.Im;
  331.     CAddF := Sum[SpSum];
  332.     IncPtr(Sum, SpSum);
  333.     end;
  334.  
  335. procedure CAdd(A, B : Complex; Result : Complex);
  336.   begin
  337.     Result^ := CAddF(A, B)^;
  338.     end;
  339.  
  340. function RxCF(A : ComplexElement; B : Complex) : Complex;
  341.   begin
  342.     CmpError := cxOK;
  343.     ChkInit;
  344.     RProd[SpRProd]^.Re := A * B^.Re;
  345.     RPRod[SpRProd]^.Im := A * B^.Im;
  346.     RxCF := RProd[SpRProd];
  347.     IncPtr(RProd, SpRProd);
  348.     end;
  349.  
  350. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  351.   begin
  352.     Result^ := RxCF(A, B)^;
  353.     end;
  354.  
  355. function CSubF(A, B : Complex) : Complex;
  356.   begin
  357.     CmpError := cxOK;
  358.     ChkInit;
  359.     Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
  360.     CSubF := Diff[SpDiff];
  361.     IncPtr(Diff, SpDiff);
  362.     end;
  363.  
  364. procedure CSub(A, B : Complex; Result : Complex);
  365.   begin
  366.     Result^ := CSubF(A, B)^;
  367.     end;
  368.  
  369. function CMulF(A, B : Complex) : Complex;
  370.   begin
  371.     CmpError := cxOK;
  372.     ChkInit;
  373.     Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
  374.     Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
  375.     CMulF := Prod[SpProd];
  376.     IncPtr(Prod, SpProd);
  377.     end;
  378.  
  379. procedure CMul(A, B : Complex; Result : Complex);
  380.   begin
  381.     Result^ := CMulF(A, B)^;
  382.     end;
  383.  
  384. function CDivF(A, B : Complex) : Complex;
  385.   var
  386.     Denom : ComplexElement;
  387.   begin
  388.     CmpError := cxOK;
  389.     Denom := B^.Re * B^.Re + B^.Im * B^.Im;
  390.     if Denom = 0.0 then begin
  391.       CmpError := cxDivByZero;
  392.       ChkInit;
  393.       Quot[SpQuot]^.Re := 0.0;
  394.       Quot[SpQuot]^.Im := 0.0;
  395.       CDivF := Quot[SpQuot];
  396.       IncPtr(Quot, SpQuot);
  397.       exit;
  398.       end;
  399.     Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
  400.     CDivF := Quot[SpQuot];
  401.     IncPtr(Quot, SpQuot);
  402.     end;
  403.  
  404. procedure CDiv(A, B : Complex; Result : Complex);
  405.   begin
  406.     Result^ := CDivF(A, B)^;
  407.     end;
  408.  
  409. procedure C2P(A : Complex; Result : Complex);
  410. {Transforms a complex in cartesian form into polar form}
  411. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  412.   begin
  413.     CmpError := cxOK;
  414.     Result^.Re := CAbsF(A);
  415.     if abs(A^.Re) > abs(A^.Im) then begin
  416.       {Use the tangent formulation}
  417.       Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
  418.       end
  419.     else begin
  420.       {Use the cotangent formulation}
  421.       Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
  422.       end;
  423.     if A^.Im > 0 then begin {first or second quadrant}
  424.       if A^.Re < 0.0 then {Second quadrant}
  425.         Result^.Im := Pi - abs(Result^.Im);
  426.       end
  427.     else begin {third or fourth quadrant}
  428.       if A^.Re > 0.0 then {Fourth quadrant}
  429.         Result^.Im := TwoPi - abs(Result^.Im)
  430.       else begin {Third Quadrant}
  431.         Result^.Im := Pi + abs(Result^.Im);
  432.         end;
  433.       end;
  434.     if Result^.Im = TwoPi then Result^.Im := 0.0;
  435.     end;
  436.  
  437. function C2PF(A : Complex) : Complex;
  438.   begin
  439.     ChkInit;
  440.     C2P(A, Polar[SpPolar]);
  441.     C2PF := Polar[SpPolar];
  442.     IncPtr(Polar, SpPolar);
  443.     end;
  444.  
  445. procedure P2C(A : Complex; Result : Complex);
  446. {Transforms a complex in polar form into cartesian form}
  447. {The magnitude is stored in A^.Re and the angle in A^.Im}
  448.   begin
  449.     CmpError := cxOK;
  450.     Result^.Re := A^.Re * cos(A^.Im);
  451.     Result^.Im := A^.Re * sin(A^.Im);
  452.     end;
  453.  
  454. function P2CF(A : Complex) : Complex;
  455.   begin
  456.     ChkInit;
  457.     P2C(A, Cartsn[SpCartsn]);
  458.     P2CF := Cartsn[SpCartsn];
  459.     IncPtr(Cartsn, SpCartsn);
  460.     end;
  461.  
  462. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  463. {Raises a complex (in polar form) to a real power, returning the result
  464.  also in polar form}
  465.   begin
  466.     CmpError := cxOK;
  467.     if A^.Re <= 0.0 then begin
  468.       CmpError := cxBadMagnitude;
  469.       ChkInit;
  470.       Power[SpPower]^.Re := 0.0;
  471.       Power[SpPower]^.Im := 0.0;
  472.       exit;
  473.       end;
  474.     Power[SpPower]^.Re := exp(R * ln(A^.Re));
  475.     Power[SpPower]^.Im := R * A^.Im;
  476.     CpPwrRF := Power[SpPower];
  477.     IncPtr(Power, SpPower);
  478.     end;
  479.  
  480. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  481.   begin
  482.     Result^ := CpPwrRF(A, R)^;
  483.     end;
  484.   end.
  485.